####### Descriptive figures ##############

# This script generates descriptive tables and figures in the order they appear in the main text and appendix.
# The generated figures are shown in the main text and in Appendix A.


# This script produces:
# Figure 2 (main text)
# Figure 3 (main text)
# Figure A2 to Figure A5 (appendix)
# Tables A2 to A4 (appendix)


# Session info
# R version 4.1.0 (2021-05-18)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS 13.1

# Load data ---------------------------------------------------------------

# define threshold of days between protest events for protest wave definition
threshold = 8

# load data
load(paste0("data/navco_week_replication",threshold,".Rda"))
load(paste0("data/navco_event_replication",threshold,".Rda"))
load("data/navco_grid_map_replication.Rda")
navco_grid_map <- navco_grid_map %>% st_as_sf()

# Mass Mobilization Data 

# load and recode data 
load("data/massmo_vdem_reign_replication.Rda") 

# data for plotting grids


# Main text tables and figures --------------------------------------------


# Figure 1 --------------------------------------------------------------

massmo_vdem_reign_summary <-  massmo_vdem_reign %>% 
  mutate(concessions_label = case_when(con_dummy == 0 ~ "No concession (n = 4,403)",
                                       con_dummy == 1 ~ "Concession (n = 479)",
                                       T ~ NA_character_)) %>% 
  group_by(concessions_label) %>% 
  summarise(mean_length = round(mean(length, na.rm = T),1)) 


# Plot densities of protest lengths (filter to keep only events longer than 1 day) by concession
# , caption = "Data source: MM (Clark/Regan 2020), sample: all events with a minimum duration of two days"
ggplot(data = massmo_vdem_reign %>% filter(dem == 0), aes(x = length, y = after_stat(density))) +
  geom_histogram(alpha=0.6, binwidth = 1, fill = "darkgreen", color = "black") +
  scale_x_continuous("Protest duration (days)", limits = c(0,30)) +
  scale_y_continuous("Share of protest intervals", limits = c(0,1), breaks = seq(0,1,.2)) +
  geom_vline(data = massmo_vdem_reign_summary, aes(xintercept=mean_length), lty = "dashed") +
  facet_wrap(~ concessions_label) +
  theme_lh +
  theme(legend.position = "none",
        strip.background = element_rect(colour="black",
                                        fill="white")) +
  geom_text(data = massmo_vdem_reign_summary, aes(x = mean_length, y = 1),
            label = paste0("avg. = ", sprintf('%.1f', massmo_vdem_reign_summary$mean_length)),
            hjust = -0.15, vjust = + 0.25)

#save as PDF
ggsave(filename = "output/figures/Figure1.pdf", width = 6, height = 3.5, dpi = "retina")


# Figure 2a ----------------------------------------------------------------

navco_event %>% 
  filter(camp_goals > 4 & camp_goals != 4) %>% 
  mutate(
    count = 1,
    no_response = ifelse(con_dummy == 0 & state_violence == 0, 1, 0)) %>% 
  group_by(camp_goals) %>% 
  summarise(count = sum(count),
            con = sum(con_dummy),
            state_violence = sum(state_violence),
            no_response = sum(no_response),
            protest = sum(protest),
            con_dummy = sum(con_dummy)) %>% 
  gather(key = "key", value = "value", con:no_response) %>% 
  mutate(camp_goals = factor(camp_goals),
         camp_goals = dplyr::recode(camp_goals, "1"="Unspecified","2"="Anti occupation", "3"="Greater autonomy", "4"="Territorial secession", "5"="Policy change", "6"="Institutional reform", "7"="Regime change"),
         percent_con = (con_dummy/protest)*100,
         key = factor(key, levels = c("no_response", "state_violence","con")),
         camp_goals = fct_reorder(camp_goals, percent_con)) %>%
  ggplot(aes(camp_goals, value, fill = key)) +
  geom_bar(position="fill", stat="identity", alpha = 0.8) +
  coord_flip() +
  labs(y = "Share of responses",
       x = "",
       fill = "") +
  scale_fill_manual(values = col3,
                    labels = c('No information','State violence','Concession')) +
  theme_lh
ggsave("output/figures/Figure2a.pdf", scale = 1, dpi = "retina", width = 18, height = 11, units = "cm")



# Figure 2b ----------------------------------------------------------------
navco_event %>% 
  filter(camp_goals > 4) %>% 
  mutate(
    count = 1,
    no_response = ifelse(con_dummy == 0 & state_violence == 0, 1, 0)) %>% 
  group_by(camp_goals) %>% 
  summarise(count = sum(count),
            con = sum(con_dummy),
            state_violence = sum(state_violence),
            no_response = sum(no_response),
            protest = sum(protest),
            con_dummy = sum(con_dummy)) %>% 
  gather(key = "key", value = "value", con:no_response) %>% 
  mutate(camp_goals = factor(camp_goals),
         camp_goals = dplyr::recode(camp_goals, "5"="Policy change", "6"="Institutional reform", "7"="Regime change"),
         
         percent_con = (con_dummy/protest)*100,
         key = factor(key, levels = c("no_response", "state_violence","con")),
         camp_goals = fct_reorder(camp_goals, percent_con)) %>%
  ggplot(aes(camp_goals, value, fill = key)) +
  geom_col(position="stack",  alpha = 0.8) +
  coord_flip() +
  labs(y = "Number of responses",
       x = "",
       fill = "") +
  scale_fill_manual(values = col3,
                    labels = c('No information','State violence','Concession')) +
  theme_lh
ggsave("output/figures/Figure2b.pdf", scale = 1, dpi = "retina", width = 18, height = 11, units = "cm")


# Figure 3 ----------------------------------------------------------------

navco_event %>% 
  group_by(move_id) %>%
  mutate(obs = row_number()) %>% 
  summarise(length = as.numeric(max(length)),
            Size = sum(size),
            con_sum = as.numeric(sum(con_dummy)),
            con_scope = max(con_scope),
            country = first(country),
            year = first(as.character(year))) %>% 
  mutate(country = capitalize(country),
         name = paste(country, year, sep = " "),
         con_scope = factor(con_scope, ordered = T),
         con_dummy = ifelse(con_sum == 0, 0, 1),
         con_dummy = as.factor(con_dummy)) %>% 
  filter(length < 50) %>% 
  ggplot(aes(x=length, con_sum, label = country, colour = con_dummy, size = Size)) +
  geom_jitter(alpha = 0.7, width = 0.2, show.legend = T) +
  geom_text_repel(aes(label=ifelse(length>40 | con_sum > 4,as.character(name),'')),hjust=0,vjust=0.5, show.legend = F, colour = "black", size = 4) +
  scale_colour_manual(values = col_con, name = "Concession",
                      labels = c("No concession", "Concession")) +
  scale_size(breaks = c(0,1000,10000, 100000, 2500000), range = c(1,10)) +
  labs(
    x = "Protest length (days)",
    y = "Sum of concessions") +
  theme_lh
ggsave("output/figures/Figure3.pdf", scale = 1, dpi = "retina", width = 26, height = 15, units = "cm")



# Appendix figures and tables --------------------------------------------------



# Figure A1 --------------------------------------------------------------


# Run negative binomial model on protest length and include protest characteristics from MM
negbin <- glm.nb(length ~ con_dummy + state_violence + protesterviolence + size_cat_f + log(protestnumber) +
                   labor_issue + land_issue + police_issue + process_issue + price_issue + politician_issue +
                   social_issue + state_violence +
                   ccode_f + year_f,
                 data = massmo_vdem_reign %>% filter(dem == 0) ,
                 control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))

# cluster standard errors
coeftest(negbin, vcov = vcovCL(negbin))

# Plot coefficients

# Note: If this plot does not reproduce, try re-installing sjPlot and its dependencies
# helpful post: https://stackoverflow.com/questions/74635326/problem-loading-sjplot-after-installing-parameters
sjPlot::plot_model(negbin, terms = c("con_dummy", "state_violence",
                             "protesterviolence", "size_cat_f1",
                             "size_cat_f2",
                             "size_cat_f3",
                             "size_cat_f4",
                             "log(protestnumber)",
                             "labor_issue1",
                             "land_issue1",
                             "police_issue1",
                             "process_issue1",
                             "price_issue1",
                             "politician_issue1",
                             "social_issue1"),
           axis.lim = c(0.1,5),
           axis.labels = c("Issue: soc. restrictions", "Issue: politician", "Issue: prices/tax policy",
                           "Issue: political process", "Issue: police brutality","Issue: land/farm",
                           "Issue: labor/wage dispute", "Protests (log)", "Size > 10000",
                           "Size 2000-10000", "Size 100-1999", "Size 50-99",
                           "Protester violence", "State violence", "Concessions"),
           vcov.args = vcovCL(negbin), 
           show.values = TRUE,
           value.offset = .35, show.p = T, value.size = 3) +
  theme_sjplot2() + 
  ggtitle("") + 
  theme(panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
  scale_color_sjplot("circus")  +
  geom_hline(yintercept = 1, color = "black", lty= "dashed")

ggsave(filename = "output/figures/FigureA1.pdf", width = 5, height = 5, dpi = "retina")


# Table A1 ----------------------------------------------------------------
negbin <- glm.nb(length ~ con_dummy + state_violence + protesterviolence + size_cat_f + log(protestnumber) +
                   labor_issue + land_issue + police_issue + process_issue + price_issue + politician_issue +
                   social_issue + state_violence +
                   ccode_f + year_f,
                 data = massmo_vdem_reign %>% filter(dem == 0) ,
                 control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))

# cluster standard errors
coeftest(negbin, vcov = vcovCL(negbin))

n_distinct(massmo_vdem_reign$ccode)
stargazer(negbin,
          type = "latex",
          float = F,
          float.env = "table",
          no.space = T,
          dep.var.caption = "Protest duration (days)",
          dep.var.labels = c(""),
          model.numbers = F,
          column.labels = c("1"),
          notes.append = T,
          notes = c("Standard errors are clustered at the level of countries."),
          covariate.labels = c("Concession (dummy)",
                               "State violence",
                               "Protester violence",
                               "Size 50-99",
                               "Size 100-1999",
                               "Size 2000-10000",
                               "Size > 10000",
                               "Protests (log)",
                               "Issue: labor/wage dispute",
                               "Issue: land/farm",
                               "Issue: police brutality",
                               "Issue: political process",
                               "Issue: prices/tax policy",
                               "Issue: politician",
                               "Issue: soc. restrictions"
          ),
          omit= c('year', 'ccode'),
          omit.stat = c("logrank", "wald", "ll", "theta"),
          add.lines = list(c("Country dummies", "Yes"),
                           c("Year dummies", "Yes"),
                           c("Number of countries", "102")),
          label = "table:massmo-negbin",
          out = paste0("output/tables/TableA1.tex"))


# Figure A2 ---------------------------------------------------------------

navco_grid_map %>%
  filter(country == "China") %>% 
  ggplot(aes(fill = I(log(events+1)))) + 
  geom_sf() + 
  scale_fill_distiller(type = "seq",
                       palette = "Blues",
                       direction = 1,
                       labels = scales::comma, 
                       name = "Events (log)") + 
  theme(panel.grid.major = element_blank()) +
  labs(title = "China")

 ggsave("output/figures/FigureA2.pdf", width = 8, height = 8, dpi = "retina")
# Table A2 ----------------------------------------------------------------

countries <- navco_week %>% 
  arrange(country, year) %>% 
  mutate(Country = capitalize(country)) %>% 
  group_by(Country) %>% 
  summarise(Protests = sum(protest),
            Concessions = sum(con_sum),
            year_first = first(year),
            year_last = last(year)
  ) %>% 
  mutate(Years = paste(year_first, year_last, sep = " - "),
         year_first = NULL,
         year_last = NULL,
         country = NULL)


print(xtable(countries,
             type = "latex",
             digits = 0,
             caption = "Case selection of protests from NAVCO 3",
             label = "tab:countries"),
             file = "output/tables/TableA2.txt")


# Figure A3a --------------------------------------------------------------

navco_event %>% 
  mutate(no_response = ifelse(con_dummy == 0 & state_violence == 0, 1, 0)) %>% 
  group_by(year) %>% 
  summarise(protest_nr = sum(protest),
            con_nr = sum(con_dummy),
            state_violence = sum(state_violence)
  ) %>%
  mutate(year = lubridate::ymd(year, truncated = 2L)) %>% 
  gather(key = "key",
         value = "value",
         protest_nr:state_violence) %>% 
  mutate(key = factor(key, levels = c("protest_nr","state_violence", "con_nr"))) %>% 
  ggplot(aes(fill=key, y=value,x = year)) +
  geom_col( position = "stack",show.legend = T,  alpha = 0.8) + 
  scale_fill_manual(values = col3,
                    labels = c("All protests", "State violence", "Concessions"),
                    name= "", position="bottom") + 
  labs(
    x="Years",
    y = '') +
  theme_lh
ggsave("output/figures/FigureA3a.pdf", scale = 1, dpi = "retina", width = 16, height = 12, units = "cm")


# Figure A3b --------------------------------------------------------------

navco_event %>% 
  mutate(no_response = ifelse(con_dummy == 0 & state_violence == 0, 1, 0)) %>% 
  group_by(year) %>%
  summarise(protest_nr = sum(protest),
            con_nr = sum(con_dummy),
            state_violence = sum(state_violence)) %>%
  gather(key = "key",
         value = "value",
         protest_nr:state_violence) %>% 
  mutate(key = factor(key, levels = c("protest_nr", "state_violence", "con_nr"))) %>% 
  group_by(year, key) %>% 
  summarise(n = sum(value)) %>% 
  mutate(percentage = n/sum(n)) %>% 
  ggplot(aes(fill=key, y=percentage,x = year)) +
  geom_col(position = "stack", alpha = 0.8) + 
  scale_fill_manual(values = col3,
                    labels = c("All protests", "State violence", "Concessions"),
                    name= "") + 
  labs(
    x="Years",
    y = '') +
  theme_lh
ggsave("output/figures/FigureA3b.pdf", scale = 1, dpi = "retina", width = 16, height = 12, units = "cm")


# Table A3 ----------------------------------------------------------------
 navco_event %>% 
  count(con_cat) 

#latex code:
# \begin{table}[H]
# \centering
# \begin{tabular}{rrrrr}
# \hline
# & Sum \\ 
# \hline
# No concession & 1678  \\ 
# Neutral &  832  \\ 
# Non-material &  25\\ 
# Material &   24  \\ 
# Full accommodation &     26  \\ 
# \hline
# \end{tabular}
# \caption{Occurrence of concession types in relation to number or protest events}
# \label{tab:con_cat}
# \end{table}

# Figure A4a --------------------------------------------------------------

navco_event %>% 
  mutate(regime = regime_reign,
         count = 1,
         no_response = ifelse(con_dummy == 0 & state_violence == 0, 1, 0)) %>% 
  group_by(regime) %>% 
  summarise(count = sum(count),
            con = sum(con_dummy),
            state_violence = sum(state_violence),
            no_response = sum(no_response),
            protest = sum(protest),
            con_dummy = sum(con_dummy)) %>% 
  gather(key = "key", value = "value", con:no_response) %>% 
  mutate(regime = factor(regime),
         percent_con = (con_dummy/protest)*100,
         key = factor(key, levels = c("no_response", "state_violence","con"))) %>%
  filter(is.na(regime) == F) %>%
  ggplot(aes(reorder(regime,percent_con), value, fill = key)) +
  geom_bar(position="fill", stat="identity", alpha = 0.8) +
  coord_flip() +
  labs(
    y = "Share of responses",
    x = "",
    fill = "") +
  scale_fill_manual(values = col3,
                    labels = c('Missing information','State violence','Concession')) +
  theme_lh
ggsave("output/figures/FigureA4a.pdf", scale = 1, dpi = "retina", width = 18, height = 12, units = "cm")


# Figure A4b --------------------------------------------------------------

navco_event %>% 
  mutate(regime = regime_reign,
         count = 1,
         no_response = ifelse(con_dummy == 0 & state_violence == 0, 1, 0)) %>% 
  group_by(regime) %>% 
  summarise(count = sum(count),
            con = sum(con_dummy),
            state_violence = sum(state_violence),
            no_response = sum(no_response),
            protest = sum(protest),
            con_dummy = sum(con_dummy)) %>% 
  gather(key = "key", value = "value", con:no_response) %>% 
  mutate(regime = factor(regime),
         percent_con = (con_dummy/protest)*100,
         key = factor(key, levels = c("no_response", "state_violence","con"))) %>%
  filter(is.na(regime) == F) %>%
  ggplot(aes(reorder(regime,percent_con), value, fill = key)) +
  geom_col(position="stack", alpha = 0.8) +
  coord_flip() +
  labs(
    y = "Number of responses",
    x = "",
    fill = "") +
  scale_fill_manual(values = col3,
                    labels = c('No information','State violence','Concession')) +
  theme_lh
ggsave("output/figures/FigureA4b.pdf", scale = 1, dpi = "retina", width = 18, height = 12, units = "cm")




# Table A4 ----------------------------------------------------------------

corr_dt <- navco_week %>% dplyr::select(
  con_sum, con_dummy, protest_end,
  # control: protest characteristic
  size_log,
  protesterviolence, goal_policy, goal_reform, goal_regime,
  # control: government characteristic
  state_violence,
  v2x_clphy, v2x_polyarchy, gdp_log, pop_log,
  #time
  diff, protests_future, protests_past
)

corr_dt <- data.frame(corr_dt)
stargazer(corr_dt,
          type = "latex",
          digits = 2,
          covariate.labels = c("Concession (sum)",
                               "Concession (dummy)",
                               "Protest end",
                               "Number of protesters (log)",
                               "Protester violence",
                               "Demand: policy change",
                               "Demand: institutional reform",
                               "Demand: regime change",
                               "State violence",
                               "Physical violence index",
                               "Democracy index",
                               "GDP pc (lag and log)",
                               "Population (lag and log)",
                               "Days since last protest",
                               "Days since last protest (squared)",
                               "Upcoming protests",
                               "Past protests"),
          float = F,
          float.env = "table",
          out = "output/tables/TableA4.tex")

# Figure A5 ---------------------------------------------------------------

supp.labs <- c("No concession", "Concession")
names(supp.labs) <- c("0", "1")

navco_week %>% 
  count(con_dummy, protests_future)%>% 
  mutate(n = n/sum(n)) %>% 
  ggplot(aes(protests_future, n)) +
  geom_col(colour = "white") +
  labs(
    x = "Protests in the next week",
    y = "Share of observations") +
  facet_wrap(~con_dummy, labeller = labeller(con_dummy = supp.labs)) +
  theme_lh
ggsave("output/figures/FigureA5.pdf", scale = 1, dpi = "retina", width = 16, height = 12, units = "cm")


